%% Manylabs robust predictive analysis
% Replicates Figure A-7
% Run "ManyLabs_main.m" first

%% Parametrization
% Number of bins for PIT histograms
bin_num = 5;

% Set max iteration number and tolerance
max_iter = 5000;
tol = 0.0001;

% Random seed
seed = 101;

% Sample size of simulated coverage & S stats
sim_cov = 5000;
sim_S = 5000;

%% Case 1: N=66, J=5 (All effects included, nu estimated by EB)
N = 66;
J = 5;
highlight_indices = [37:42, 49:60];

% Extract study estimates and ses
sheet1 = 'N=66, J=5';
estimates1 = xlsread('../Data/ManyLabs_Data.xlsx', sheet1, 'C3:H68'); 
se1 = xlsread('../Data/ManyLabs_Data.xlsx', sheet1, 'N3:S68');   

% Randomly assign baseline and validation studies
[hat_theta_1, hat_vartheta_1, rep_base_se2_1, rep_val_se2_1] = random_assignment(estimates1, se1, seed);

% Estimate nu from pooled data
EB_est_nu_1 = MLE_nu(hat_theta_1, rep_base_se2_1, 0, max_iter, tol);

% Check whether the validation study estimates are located in the predicted interval
cov_stat_sim = simulate_cov_stat_distribution(N, sim_cov);
[bar_tau_1, predict_se_1, empirical_cov_1, cov_p_value_1, cov_abs_dev_1] = check_CI_robust(hat_theta_1, hat_vartheta_1, rep_base_se2_1, rep_val_se2_1, EB_est_nu_1, cov_stat_sim);

xtick = 0:10:70;
ytick = -4:1:3;
interval_1 = plot_intervals_fixord(hat_vartheta_1, bar_tau_1, predict_se_1, highlight_indices, order_1, xtick, ytick);
saveas(interval_1, "../Results/Robust_EB_Interval_N_66_J_5.png")


%% Case 2: N=132, J=2 (All effects included, nu estimated by EB)
N = 132;
J = 2;
highlight_indices = [73:84, 97:120];

% Extract study estimates and ses
sheet2 = 'N=132, J=2';
estimates2 = xlsread('../Data/ManyLabs_Data.xlsx', sheet2, 'C3:E134'); 
se2 = xlsread('../Data/ManyLabs_Data.xlsx', sheet2, 'K3:M134');  

% Randomly assign baseline and validation studies
[hat_theta_2, hat_vartheta_2, rep_base_se2_2, rep_val_se2_2] = random_assignment(estimates2, se2, seed);

% Estimate nu from pooled data
EB_est_nu_2 = MLE_nu(hat_theta_2, rep_base_se2_2, 0, max_iter, tol);

% Check whether the validation study estimates are located in the predicted interval
cov_stat_sim = simulate_cov_stat_distribution(N, sim_cov);
[bar_tau_2, predict_se_2, empirical_cov_2, cov_p_value_2 , cov_abs_dev_2] = check_CI_robust(hat_theta_2, hat_vartheta_2, rep_base_se2_2, rep_val_se2_2, EB_est_nu_2, cov_stat_sim);

xtick = 0:20:140;
ytick = -4:1:3;
interval_2 = plot_intervals_fixord(hat_vartheta_2, bar_tau_2, predict_se_2, highlight_indices, order_2, xtick, ytick);
saveas(interval_2, "../Results/Robust_EB_Interval_N_132_J_2.png")


%% Case 3: N=48, J=5 (three irreplicable studies excluded, nu estimated by EB)
N = 48;
J = 5;
highlight_indices = [];

% Extract study estimates and ses
sheet3 = 'N=48, J=5';
estimates3 = xlsread('../Data/ManyLabs_Data.xlsx', sheet3, 'C3:H50'); 
se3 = xlsread('../Data/ManyLabs_Data.xlsx', sheet3, 'N3:S50');   

% Randomly assign baseline and validation studies
[hat_theta_3, hat_vartheta_3, rep_base_se2_3, rep_val_se2_3] = random_assignment(estimates3, se3, seed);

% Estimate nu from pooled data
EB_est_nu_3 = MLE_nu(hat_theta_3, rep_base_se2_3, 0, max_iter, tol);

% Check whether the validation study estimates are located in the predicted interval
cov_stat_sim = simulate_cov_stat_distribution(N, sim_cov);
[bar_tau_3, predict_se_3, empirical_cov_3, cov_p_value_3, cov_abs_dev_3] = check_CI_robust(hat_theta_3, hat_vartheta_3, rep_base_se2_3, rep_val_se2_3, EB_est_nu_3, cov_stat_sim);

xtick = 0:10:50;
ytick = -4:1:3;
interval_3 = plot_intervals_fixord(hat_vartheta_3, bar_tau_3, predict_se_3, highlight_indices, order_3, xtick, ytick);
saveas(interval_3, "../Results/Robust_EB_Interval_N_48_J_5.png")


%% Case 4: N=96, J=2 (three irreplicable studies excluded, nu estimated by EB)
N = 96;
J = 2;
highlight_indices = [];

% Extract study estimates and ses
sheet4 = 'N=96, J=2';
estimates4 = xlsread('../Data/ManyLabs_Data.xlsx', sheet4, 'C3:E98'); 
se4 = xlsread('../Data/ManyLabs_Data.xlsx', sheet4, 'K3:M98');   

% Randomly assign baseline and validation studies
[hat_theta_4, hat_vartheta_4, rep_base_se2_4, rep_val_se2_4] = random_assignment(estimates4, se4, seed);

% Estimate nu from pooled data
EB_est_nu_4 = MLE_nu(hat_theta_4, rep_base_se2_4, 0, max_iter, tol);

% Check whether the validation study estimates are located in the predicted interval
cov_stat_sim = simulate_cov_stat_distribution(N, sim_cov);
[bar_tau_4, predict_se_4, empirical_cov_4, cov_p_value_4, cov_abs_dev_4] = check_CI_robust(hat_theta_4, hat_vartheta_4, rep_base_se2_4, rep_val_se2_4, EB_est_nu_4, cov_stat_sim);

xtick = 0:20:100;
ytick = -4:1:3;
interval_4 = plot_intervals_fixord(hat_vartheta_4, bar_tau_4, predict_se_4, highlight_indices, order_4, xtick, ytick);
saveas(interval_4, "../Results/Robust_EB_Interval_N_96_J_2.png")


%% Summarize the EB results
Robust_EB_empirical_cov = [empirical_cov_1, empirical_cov_2, empirical_cov_3, empirical_cov_4]
Robust_EB_cov_p_value = [cov_p_value_1, cov_p_value_2, cov_p_value_3, cov_p_value_4]
Robust_EB_est_nu = [EB_est_nu_1, EB_est_nu_2, EB_est_nu_3, EB_est_nu_4]
